Comparison of two different implementations of a 
finite-difference-method for first-order pde in 
Mathematica® and Matlab® 
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| 7 | \ Abstract 

In this article two implementations of a symmetric finite difference algorithm (ftcs-method) 
for a first-order partial differential equation are discussed. The considered partial differen- 
tial equation discribes the time evolution of the crack length ditribution of microcracks in 
brittle materia. 
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1 Physics and analytical solution 



The growth rate of microcracks in a brittle material can be discribed by a mesoscopic equation, 
i/-) \ Here the specialized version for uniaxial loading is presented. 

8; df{i,t) _ i9/ 2 ^(/,t)/(/,_t) 



dt I 2 dl 

f(l, t) is the distribution function for the crack length I at time t, V\ — I is the growth velocity 
q ! of the cracks. A Rice-Griffith-like dynamic is assumed for crack growth, which gives 

/ -a' + p'laitf , if a' < /37a 2 
, I — s (2) 

\_j . I , otherwise 

The theory is given in detail in CJ. For an exponential (or a step- wise) initial condition and 
constant loading speed it is possible to give an exact analytical solution, which is also presented 
in HI and looks like: 



f(l,t) 
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f(l,0) , otherwise. 
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(a) step-wise initial condition (b) step-wise initial condition 




(c) exponential initial condition 
Figure 1: Analytical solution for crack lenght distributuion. Shown is f(l, t)l 2 over / and t. 
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2 The numerical algorithm 

The partial differential equations are first order in time and crack length. Two different al- 
gorithms, upwind and fcts (forward time centered space), have been tested. The symmetric 
algorithm (ftcs-method) is in this case a bit more stable than the upwind, which is somewhat 
astonishing. Both algorithms can be found in 0. 

In the symmetric algorithm /(/, t + 1) is calculated from f(l — 1, t), /(/, t) and f(l + 1, t) by 



This is what the main part of the implementation in Mathematica® looks like: 

For[t = 1, t < TMAX, t++, For[l = 1, 1 < LMAX, 1++, 
If [l*dl* (t*dt) A 2*be - al > 0, 
If[l == 0, f[l, t + 1] = f[l, t], 
f[l, t + 1] = 

f[l, t]*(l - (3*be* (t*dt) A 2 - (2*al) / (dl*l) ) *dt) - 
(dl*l*be* (t*dt) A 2 - al) *dt/dl* (f [1 + 1, t] - f[l - 1, t] ) ] , 
f[l, t + 1] = f[l, t] 
] 

] 



This is what the main part of the implementation in Matlab® looks like: 

for t=l:l:TMAX-l 
for 1=1:1:LMAX-1 

if (1 . *dl . *sigma (t) . *beta-alpha > 0) 

f (1, t + 1) =f (1, t) . * (1- (3. *beta. * (sigma (t) ) . A 2 - ... 
(2 . *alpha) . / (dl. *1) ) . *dt) - (dl. *1. *beta. * (sigma (t) ) . A 2 - ... 
alpha) . *dt. / (2*dl) . * (f (1 + 1, t) -f (1-1, t) ) ; 
else 

f (l,t + l)=f (l,t) ; 
end 




x— (f(l + l, t )-f(l-l,t)) 



(4) 



end 



end 
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3 Results obtained by Mathematica® 

Mathematica® is a general purpose computer algebra system (cas) by Wolfram Research 
Inc.. As one can see the solution shows some wave-like efects, which can be interpreted as a 




(a) step-wise initial condition (b) exponential initial condition 

Figure 2: Numerical solution for crack lenght distributuion calculated with Mathematica' 
sign for numerical instability. This instability is a result of the discontinous initial condition. 
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4 Results obtained by Matlab® 

Matlab® is a tool for numerical mathematics, especially designed for matrix manipulation, 
by The Math Works Inc.. Here in both cases huge instabilities occur. As soon as some cracks 




(a) step- wise initial condition 



(b) exponential initial condition 



Figure 3: Numerical solution for crack lenght distributuion calculated with Matlab® 

start growing the numerical error goes to infinity. The following pictures show an extract of the 
above pictures. 





(a) step-wise initial condition 



(b) exponential initial condition 



Figure 4: Numerical solution for crack lenght distributuion calculated with Matlab® (extract 
of figure El) 
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5 Comparison of the results and conclusion 



Obviously the results, obtained with both programs, show huge errors. But in contrast to the 
results calculated with Matlab® , one can use the results obtained with Mathematica® at 
least for some rough predictions. 

The difference in the two software packages, used for these simulations, is that Matlab® uses 
floating-point variables of precision "double" (16 byte), whereas Mathematica® is capable 
of both numerical and symolic computation. Therefore it is possible that Mathematica® uses 
a much higher precision to perform some of the operations than MATLAB® . 
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